package NetWork;

import odeToJava.modules.ODE;

/**
 *
 * @author citizen5589 & MawakMan
 */
public class Omega implements ODE {
    // Iapp time var
    public static int temp = 0;
    
    // Time delay  for neurons
    double epsilon = 0.04D;
    
//    double epsilon2 = 0.04D;
//    double epsilon3 = 0.04D;
//    double epsilon4 = 0.04D;
    
     // Constant Voltage for neurons
    double I1 = 1.05D;
    double I2 = 1.02D;
    double I3 = 1.02D;
    double I4 = 1.02D;
    
    // Time separation parameter for Ca2+ oscillation
    double epsilonc = 0.04D;
   
    // Time delay for synapse ??? to be fixed
    double tau_s = 10.0D;
    // Characteristic time for Ca2+ oscillation
    double tau_c = 8.0D;
    // Time Control parameter for mediator production (Sm)
    double tau_Sm = 100.0D;
    // Time Control parameter for mediator production (Gm)
    double tau_Gm = 5.0D * tau_s;

    // Constants Values for Z Activation & relaxation
    double Ss = 1.0D;
    double hs = 0.0D;
    double ds = 3.0D;
    
    // Controls parameters for mediator production (Sm)
    double SSm = 100.0D;
    double hSm = 0.45D;
    double dSm = 3.0D;
    // Controls parameters for mediator production (Gm)
    double SGm = 100.0D;
    double hGm = 0.5D;
    double dGm = 3.0D;
    // Initial value of Z computed at program launch
    double z0 = (2.0D * ds) / (1 + 2.0D * ds + Math.exp(2.0D * I1));
    // initial state of calcium oscillator
    double r = 0.31D;
   
    // Constants values for Calcium dynamic of astrocyte 
    double c1 = 0.13D;
    double c2 = 0.9D;
    double c3 = 0.004D;
    double c4 = 2.0D / epsilonc;
    
    // postsyn neuron to astrocyte influence
    public static double alpha;

    // syn to astrocyte influence
    public static double beta;
    
    // AMPA from astrocyte
    private static double gamma;
    
     // GABA from astrocyte
    private static double delta;
    
    // Conductivity parameter for Isyns
    // n1 synapses (0.2)
    double ks12 = 0.2D;
    double ks13 = 0.2D;
    double ks14 = 0.0D;
    //double ks15 = 0.2D;
    // n2 synapses 
    double ks23 = 0.2D;
    //double ks24 = 0.2D;
    //double ks25 = 0.2D;
    // n3 synapses
    double ks34 = 0.0D;
    //double ks35 = 0.2D;
    // n4 synapses
    double ks45 = 0.2D;
    
    // Applied Voltage on presyn neuron
    double Iapp = 0.0D;
    
    private static double s = Simulation.getStep();
    private static double timeFactor = 4/s;
    private static double stimStart;
    private static double stimEnd;
    private static double stimStart0;
    private static double stimEnd0;
    
    
    public double[] f(double t, double[] x) {

        // Apply Iapp
        temp++; 
        if ( ((temp>stimStart) && (temp<stimEnd)) || ((temp>stimStart0) && (temp<stimEnd0)) )
        {Iapp = 0.8D;} else {Iapp=0.0D;}
        
        double[] dx = new double[x.length];
        
        // astrocyte factor in synapses
        double deltaGm = delta * x[17];
        // Iglion
        double Iglion = gamma * x[17];

        // v1
        dx[0] = (x[0] - Math.pow(x[0], 3) / 3.0D - x[1]) / epsilon;
        // w1
        dx[1] = x[0] + I1 - Iapp;
        // z1
        dx[2] = ((1.0D + Math.tanh(Ss * (x[0] - hs))) * (1.0D - x[2]) - (x[2] / ds)) / tau_s;
        // Isin1
        double Isin12 = (ks12 - deltaGm)*(x[2]-z0);
        
        // v2
        dx[3] = (x[3] - Math.pow(x[3], 3) / 3.0D - x[4]) / epsilon;
        // w2
        dx[4] = x[3] + I2 - Isin12 - Iglion;
        // z2
        dx[5] = ((1.0D + Math.tanh(Ss * (x[3] - hs))) * (1.0D - x[5]) - (x[5] / ds)) / tau_s;
        
        // v3
        dx[6] = (x[6] - Math.pow(x[6], 3) / 3.0D - x[7]) / epsilon;
        // w3
        dx[7] = x[6] + I3 - Iglion;
        //dx[7] = x[6] + I3 - Isin13 - Isin23  - Iglion;
        // z3
        dx[8] = ((1.0D + Math.tanh(Ss * (x[6] - hs))) * (1.0D - x[8]) - (x[8] / ds)) / tau_s;   
        // Isin3
        double Isin34 = (ks34 - deltaGm)*(x[8]-z0);
        
        // v4
        dx[9] = (x[9] - Math.pow(x[9], 3) / 3.0D - x[10]) / epsilon;
        // w4
        //dx[10] = x[9] + I3 - Isin14 - Isin24 - Isin34 - Iglion;
        dx[10] = x[9] + I3 - Isin34 - Iglion;
        //dx[10] = x[9] + I3 - Isin14 - Isin34 - Iglion;
        // z4
        dx[11] = ((1.0D + Math.tanh(Ss * (x[9] - hs))) * (1.0D - x[11]) - (x[11] / ds)) / tau_s;   

        // c : 12
        dx[12] = (-x[12] - (c4 * newf(x[12], x[13])) + (r + (alpha*x[4]) + (beta*x[14])
                                                          + (alpha*x[7]) + (beta*x[15])
                                                          + (alpha*x[10])+ (beta*x[16]))) / tau_c;
        // ce : 13
        dx[13] = newf(x[12], x[13]) / (epsilonc*tau_c);
        // Sm (with syn entries ) : 14,15,16
        dx[14] = ((1.0D + Math.tanh(SSm * ( x[2] - hSm))) * (1.0D - x[14]) - (x[14] / dSm)) / tau_Sm;
        dx[15] = ((1.0D + Math.tanh(SSm * ( x[5] - hSm))) * (1.0D - x[15]) - (x[15] / dSm)) / tau_Sm;
        dx[16] = ((1.0D + Math.tanh(SSm * ( x[8] - hSm))) * (1.0D - x[16]) - (x[16] / dSm)) / tau_Sm;
        // Gm : 17
        dx[17] = ((1.0D + Math.tanh(SGm * (x[12] - hGm))) * (1.0D - x[17]) - (x[17] / dGm)) / tau_Gm;
        return dx;
    }

    public double[] g(double arg0, double[] arg1) {
        double[] event = new double[1];

        return (event);
    }

    private double newf(double c, double ce) {
        double c_2 = Math.pow(c, 2);
        double ce_2 = Math.pow(ce, 2);
        double c2_4 = Math.pow(c2, 4);
        double c_4 = Math.pow(c, 4);
        return ((c1 * (c_2 / (1.0D + c_2))) - (ce_2 / (1.0D + ce_2)) * (c_4 / (c2_4 + c_4)) - (c3 * ce));
    }
    
    public static void initValues(double gammaV,double deltaV,double alphaV,double betaV,double stimStartV,double stimEndV,double stimStart0V,double stimEnd0V){
        gamma = gammaV;
        delta = deltaV;
        alpha = alphaV;
        beta = betaV;
        stimStart = stimStartV*timeFactor;
        stimEnd = stimEndV*timeFactor;
        stimStart0 = stimStart0V*timeFactor;
        stimEnd0 = stimEnd0V*timeFactor;  
        temp = 0;
    }
}
